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ABSTRACT 

There is mounting evidence for the coevolution of galaxies and their embedded massive black holes (MBHs) 
in a hierarchical structure formation paradigm. To tackle the nonlinear processes of galaxy - MBH interaction, 
we describe a self-consistent numerical framework which incorporates both galaxies and MBHs. The high- 
resolution adaptive mesh refinement (AMR) code Enzo is modified to model the formation and feedback of 
molecular clouds at their characteristic scale of 15.2 pc and the accretion of gas onto a MBH. Two major 
channels of MBH feedback, radiative feedback (X-ray photons followed through full 3D adaptive ray tracing) 
and mechanical feedback (bipolar jets resolved in high-resolution AMR), are employed. We investigate the 
coevolution of a 9.2 x 10 11 M Q galactic halo and its 1O 5 M embedded MBH at redshift 3 in a cosmological 
ACDM simulation. The MBH feedback heats the surrounding ISM up to 10 6 K through photoionization and 
Compton heating and locally suppresses star formation in the galactic inner core. The feedback considerably 
changes the stellar distribution there. This new channel of feedback from a slowly growing MBH is particularly 
interesting because it is only locally dominant, and does not require the heating of gas globally on the disk. 
The MBH also self-regulates its growth by keeping the surrounding ISM hot for an extended period of time. 
Subject headings: galaxies: formation — stars: formation — galaxies: active — galaxies: nuclei 



1. INTRODUCTION 

Ever since the discovery of the ubiquitous existence of su- 
permassive black holes at the c enters of massive galaxies (e.g. 
iKormendv & Richstond [T995). a plethora of evidence has ac- 
cumulated to indicate the coevolution of galaxies and their 
embedded massive black holes (MBHs). The observed tight 
correl ation between MBH masses and bulge velocity disper- 
sions dFerrarese & Merrittl 120001 iGebhardt et al.l 2000) have 
bolstered the idea that the fates of a host galaxy and its embed- 
ded MBH are fundamental ly intertwined and h eavily affected 
by each other's influence (Silk & Rees 1998; Kauffmann & 
Haehnelt 2000 UWyithe & Loebll2003rK 

Recent observations provide more solid constraints on the 
coevolution of galaxies and MBHs. For example, cosmologi- 
cal star formation history and black hole accretio n history are 
measu red to be proportional to each other (e.g. IZheng et al.1 
2009). Me rging of galaxies is b elieved to induce quasar ac- 
tivity (e.g. Hopkins et al. 2008), and the existence of high- 
redshift quasars dFan et al.ll2006f) indicate the rapid growth of 
black hole masses in the early phase of hierarchical structur e 
formation, most likely by mergers (Haima n~& Loebll2001l) . 
Unmistakably it is a complicated and highly nonlinear pro- 
cess for a galaxy to affect its embedded MBH, and vice versa. 
Therefore, developing a numerical tool which incorporates 
both galaxies and MBHs in one self-consistent framework is 
indispensible to fully co mprehend thei r coevo lution. 

The seminal work by Spring eTet al.l (l2005bh to include ac- 
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cretion and feedback of a MBH in a galactic simulation has 
been followed by many detailed investigations. These stud- 
ies have helped extend our understanding of galaxy - MBH 
interaction in various contexts and scales: (a) Merging of 
Milky Way sized galaxies was simulated to show that quasar- 
like MBH feedback drives a massive gas outflow leading to 
quenche d star formation, and to the observed Mr h — Cbuige 
relation dSpringel et alj|2005allbl: iDi Matteo et al.[|2005b Jo- 
hansson et al. l2009h . (b) Successive mergers of galaxies and 
MBHs were performed in a cosmological volume to yie ld a 
viable route to form high-redshift quasars (Li et al. 2007; Si- 
jacki et al. 2009|). (c) MBH feedback at the center of a galaxy 
cluster was demonstrated to rel ease sufficient energ y to stop 
an overlv cooled inflow of gas (Siiacki et al. 2007; Booth & 
Schave T2009l;lTevssier et al.ll2010l : lDubois et al.ll2010l) . 

Nonetheless, a comprehensive numerical understanding 
which incorporates both galaxies and MBHs is still miss- 
ing, for various reasons. First and foremost, simulated galax- 
ies do not match some of the most obvious aspects of ob- 
served galaxies. For example, simu lated galaxies are prone to 
lock baryons into too many stars dGuo et alJl2010L and ref- 
erences therein), or contain bulge-dominated disks that are 
too centrally concentrated and have a gr eatly reduced angular 
momentum relative to those observed (van den Boschll200ll : 
iKaufmann et all l2007t iPiontek & Steinmetzll2009ab . These 
problems are somewhat alleviated by loweri ng star formation 
efficiency and/or increasing stellar feedback ( Governat o et al.1 
l2007LlPiontek & Steinmetz 2009b; Agertz et alJ201 lh . or even 
by introducing a new powerful energy source such as MBH 
feedback. However, the former fix has not been entirely suc- 
cessful even with varied feedback parameters while the latter 
almost always powers large-scale gas outflow leaving behind 
a "red and dead" galaxy devoid of gas for a long time (Bor- 
gani et al.l2 004; Kra vtsov et al.ll2005HSpringel et alJ l2005b: 
iTeyssier et alJl201QT) . Obviously numerical simulations are 
still missing one or more essential ingredients. It could be the 
ignored physical processes such as stellar UV radiation and 
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magnetic fields. Or it could be the inaccurate descriptions of 
MBH accretion and feedback. 

Second, most numerical studies to date lack necessary res- 
olution and technique to describe how gas falls onto a central 
MBH and how the energy input of MBH feedback is deposited 
to its surrounding gas. While the 1 - 100 kpc resolution in 
large-scale simulations is clearly insufficient to adequately de- 
scribe the accretion flow onto a MBH, even galactic scale sim- 
ulations do not generally resolve the Bondi radius (See 32.61 
Eq.dlOb). which is required in order to trace how a MBH grav- 
itationally influences its surroundings and how the radiation 
and outflows from the MBH are thermally coupled to the gas. 
Indeed, poor resolution has forced simulators to skip the ther- 
malization process below the resolution limit, and to simply 
thermodynamically deposit MBH feedback energy near the 
MBH. While crude, it has been an effective approximation 
characterizing MBH feedback on a resolved scale (Springel 
et al. l2005bh . And it might be a fairly reasonable choice if 
MBH feedback is powerful enough to drive thermal shock 
waves (so-called "quasar-mode"; Mbh > 0.02MEdd)- How- 
ever, it can not adequately describe the energy coupling of 
the radiation from a weak, quiescent MBH ("radio-mode"; 
Mbh < 0.02M Edd ; ICroton et alj|2006t iMcNamara & Nulsenl 
120071) . For this reason injecting thermal energy in a small vol- 
ume of poorly resolved interstellar medium (ISM) can hardly 
be an accurate description of MBH feedback (See 32.6| for de- 
tailed discussion). Modeling how MBH feedback energy is 
actually coupled to the gas is a critical missing piece in con- 
temporary galaxy formation simulations. 

Third, partly due to the lack of proper resolution, most nu- 
merical calculations to date have modeled stars and MBHs 
with phenomenologically parametrized ad hoc formulations. 
Most notably, the Eddington-limited Bondi- Hoyle accretion 
estimat e employed by many authors (e.g. ISpringel et alj 
l2005bl;lDi Matteo et al.ll2005[|Johansson et al.ll2009l See 3231 
for definitions of variables) has had to be empirically boosted 
by an efficiency parameter a =10 - 300. 

M BH = m inf 47raG2 fBHPB/^BH^ P y 
V c s ^°tC ) 

While this nondimensional boost factor a is to correct the 
large-scale averaged, and probably underestimated p B near 
the MBH, a is typically fixed after the MBH has grown 
so the Bondi radius is resolved even with coarse resolu- 
tion. 8 Another example of introducing tunable parameters 
based on unknown physics is to use two different implemen- 
tations of MBH feedback, depending on the estimated accre- 
tion rate: quasar-mode feedback and rad io-mode feedback 
dSiiacki et al.ll2007t|Puchwein et al.ll2008l) . While useful in 
some applications, these ad hoc approaches ironically demon- 
strate that the physics of MBHs has not yet been adequately 
described in simulations. 

In order to circumvent the limitations of previous ap- 
proaches outlined above and to follow the actual physical 
processes between gas, stars, and MBHs, we develop a fully 
self-consistent galaxy formation simulation integrating the 
growths of both galaxies and MBHs in one comprehensive 
framework. We limit the use of ad hoc formulation but instead 
more accurately model the physics in all aspects of galaxy 
formation, namely: (a) molecular cloud formation, (b) stellar 

8 For more discussion on the boost factor a, see Johansson et al. (2009) 
or Booth & Sc have (2009). For a discussion on the accretion rate adopted in 
this work, see !]2.6i and Appendix lAl 



feedback, (c) MBH accretion, and (d) MBH feedback. Our 
code models the formation and feedback of m ol ecula r clouds 
at their characteristic scale of 15.2 pc ( 32.41 to l23T l and the 
accretion of gas onto a MBH ( 32.61 l. Two major channels 
of MBH feedback are also considered: radiative feedback 
(monochromatic X-ray photons f ollow ed through full three 
dimensional adaptive ray tracing; H2.1\ and mechanical feed- 
back (bipolar jets resolved in high-resolution adaptive mesh; 
32.8l i. We then investigate the evolution of a 9.2 x 10 n M® 
galactic halo with an embedded seed MBH of 10 5 M Q at z ~ 3 
in a cosmological ACDM simulation. 

This paper will be the first in a series that assembles a num- 
ber of high-resolution galaxy formation simulations with self- 
consistently modeled stars and MBHs. This article is orga- 
nized as follows. The physics of galaxy formation in our code 
is the topic of Ej2j followed by the initial condition of our sim- 
ulation in ij3] Slis devoted to the results of our experiments, 
with an emphasis on the feedback-regulated star formation 
and black hole growth. Discussed in 33] are the conclusions 
and the future work. 

2. MODELING THE PHYSICS OF GALAXY 
FORMATION 

The high-resolution Eulerian adaptive mesh refinement 
(AMR) code Enzo-2.0 (http://enzo.googlecode.com/; Bryan 
& Norman 1 1 997t iNorman & Brvan[l l999l;lBrvan et alj |2001l; 
lO'Shea et all 12004 INorman et al.N2007h captures the grav- 
itational collaps e of turbulent fragmentation with hig h spa- 
tial resolution dWise & Abell 1200% IWise et al.l 120081: Turk 
et al. l2009h and attains multiphase gas dynamics i n the ISM 
as it sharply resolves shocks and phase bound aries ( Sl vz et all 
12001 lAgertzltli]l2007l; iTasker et al.ll2008l) . Our enhanced 
version of Enzo contains all relevant f eatures previously dis- 
cussed in simulating galaxy evolution (|Tasker & Brvanll2006l 
120081: iKim et al.ll2008ll2009h as well as a treatment of several 
new physical processes discussed in detail in the following 
sections. 

2.1. Hydrodynamics and Gravitational Dynamics 

The ZEUS astrophysical hydrodynamics module included 
in Enzo is employed to solve the Euler equa tions for the col- 
lisional baryon fluid represented by grids dStone & Normanl 
ll992allb1:lAnninos & Normanlll994h . While known to intro- 
duce spurious effects, this scheme is widely used with AMR 
because of the stability of its solutions, and the acceptable er- 
ror when combined with high resolution. 

Dark matter, stars, and MBHs are treated as collisionless 
particles which interact only by the gravitational force. To 
evolve the particle positions and velocities, the gravitational 
dynamics are solved by an N-body adaptive particle-mesh 
solver. After particles are gridded onto the mesh by the cloud- 
in-cell interpolation, the Poisson equation is solved on the dis- 
cretized density grids via fast Fourier transform and mu ltigrid 
solvers dHocknev & Eastwoodll988tlO'Shea et aill2004h . 

2.2. Refinement Strategy 

Enzo decides whether each parent cell needs to be refined 
into eight child cells based on the mass of the cell in gas or in 
particles. The timestep is also adaptively determined level by 
level so that the timestep dt satisfies 

rfK0.3x^ = 0.3 x (sound crossing time) (2) 

for all the cells at that level. Here c s is the sound speed of 
the gas, and we choose the Courant-Friedrichs-Lewy (CFL) 
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safety number of 0.3. As decisions for refinement are made 
recursively, the resulting dataset is a nested grid-patch struc- 
ture. In our work, the grids are adaptive ly refined down to 
15.2 pc resolution. This value is in accord with the Jeans 
length for a dense gas clump of n = 125 cm~ 3 at —200 K, 
at which point a corresponding Jeans mass of 16000M© col- 
lapses to spawn a molecular cloud particle. 

We refine the cells by factors of 2 in each axis, on gas and 
particle overdensities of 8. The mass thresholds, M re f, above 
which a cell refines are functions of a refinement level / as 

<f, gas =2- 3 ™M« f , gas = 2-°- 37 « .0.125i2 bP0 A* 3 (3) 
Alport =2- al05/ < f ,p art = 2-°- 105/ -0.125i2 m p Ax 3 (4) 

where factors 0.125 = 8(l/2 3 ) 2 guarantees to refine all the 
cells of the first two nested levels f fl3.ll . Ax is the cell size 
at a root grid, and pQ = 3H 2 /%nG is the critical density. 
£2 m = 0.27, £2 b = 0.044, and H = 7 1 km s~ 'Mpc~ 1 are matter 
density, baryon density, and the Hubble constant, respectively. 

For example, at the finest static level, I = 2, a cell is refined 
if it has more mass than 8.9 x 10 5 M Q in gas or 6.7 x 10 6 M© 
in particles. At level I = 1 1 (Ax = 15.2 pc at z = 3) a cell is re- 
fined if more than 8.4 x 10 4 M@ = 5 M Je ans(125 crrT 3 ,200 K) 
in gas or 3.5 x 10 M® = 47 MoM.smaiiest in particles. This 
way we refine the grids more on small scales, which allows 
us to focus our computational resources more on the dense 
star forming regions, mak ing the simulation super-Lagrangian 
JO' Shea & Normaril2008l) . 

2.3. Chemistry and Radiative Cooling 

We use non-equilibrium chemistry model to track six 
species (H, H + , He, He + , He ++ , e~) by following six 
collisional processes among them. At the same time, 
Enzo's cooling module considers collisional excitation cool- 
ing, collisional ionization cooling, recombination cooling, 
Bremsstrahlung cooling, and CMB Compton c ooling to com- 
pute the radiative loss o f internal gas energy (lAnninos et al.l 
Il997t lAbel et al.lll997h . Added to these primordial cooling 
rates is the metallicity-dependent metal cooling rate AA(Z) = 
A n et(Z) — A net ( 0) above 10 4 K, where A net i s the net cooling 
rate tabulated in lSutherland & Dopital Jl993h . Cooling below 
10 4 K is also e nabled with fine structure metal-line coolin g 
by C, O, and Si JGlover & JappserJ2007tlWise & Abelll200l . 
This treatment ensures that a thin galactic disk forms by being 
cooled below 10 4 K, the approximate virial temperature of the 
ISM in a galactic disk. 

We further refine our module with photoionization heating 
at z < 3 by the metagalactic backgroun d UV of quasars and 
galaxies (lHaardt & Madaulll996l 12001). which is known to 
give rise to a war m diffuse ISM and preven t star formation in 
optically thin gas JCeverino & K lypin 2009). An approximate 
self-shielding fac tor is applied when the heating term is added 
JCen et al.ll2005l) . While not introducing a marked difference 
in overall results analyzed here, inclusion of this additional 
heating term results in a more realistic interstellar medium. 

2.4. Molecular Cloud Formation 

Our molecular cloud particle formation is based on Cen & 
Ostriker (1 19921) formalism with several important modifica- 
tions. With a fixed formation efficiency of e* = 0.5, the finest 
cell of physical size Ax = 15.2 pc and gas density p gas pro- 
duces a molecular cloud particle of mass 

M M c = £*P gas Ar 3 (5) 




FIG. 1. — A schematic view of the molecular cloud formation and the 
stellar feedback. 12?^,, after a molecular cloud particle (Mmc — 8000 M .-,) is 
formed, only 20% of its mass remains as an actual stellar mass M stllI (t) while 
the rest 80% has returned to the gas along with thermal feedback energy. 



when the following four criteria are met: 

(a) the proton number density exceeds the threshold 

«thres = 125 cirr 3 , 

(b) the velocity flow is converging; i.e. V • v < 0, 

(c) the cooling time f coo i is shorter than the dynamical time 
f dyn of the cell: E int /E < [37r/(32Gp gas )] 1/2 , and 

(d) the particle produced has at least M t h res = 8000 M®. 

The consequence of our criteria is the following. The gas in 
the finest cell is converted into a particle as soon as the cell 
has accumulated more than M t h res /e* = 16000 M®, the Jeans 
mass at n = 125 cirr 3 at -200 K. Because 8000 - 42000 M® 
is instantly removed from the cell every time a particle is cre- 
ated, the gas mass in the finest cell never reaches the refine- 
ment threshold Mjj 1 * = 84000 M® described in fl2~2l ensur- 
ing the consistency between the refinement criteria and the 
particle formation. 9 The values used here are in good agree- 
ment with those corresponding to collaps ing Giant Molecular 
Clouds (GMC: lMcKee & Ostrikerll2007l) where star-forming 
molecular clumps are enshrouded by cold atomic gas. 

As an additional note, differences f rom more traditional star 
particle formation criteria such as in Task er & Brvar] J2008) 
include: (a) the Jeans condition p gas Ax 3 > Mj ean s is removed 
because this condition could have allowed mass greater than 
Cleans t° accumulate while not being properly resolved un- 
til a particle finally for ms, (b) the factor Af/f<jyn in Eq.(l) of 
Taske r~& Bryanl J2008I) is removed in order to instantly create 
a particle and not leave any unresolved mass behind, and (c) 
stochastic star formation is not imposed. 

With these modifications, our criteria guarantee that a par- 
ticle forms before an unphysically large mass begins to ac- 
crete onto any unresolved dense clump. It is worth to em- 
phasize the differences between our molecular cloud forma- 
tion criteria and prior studies. While many previous studies 
with particle-based codes (e.g.lMihos et al.Hl 99 lhlKatzll 1992k 
iMihos & Hernquist|[T994l: ISpringel & Hernquistll2003f: Gov- 
ernato et al. 120071) place a star particle using the Schmidt 
relation (Psfr ~ Pgas^ ISchmidtl 1 19591) . we deposit a parti- 
cle when a gas cell of a typical molecular cloud size actu- 
ally becomes Jeans unstable. For this reason, the particle 

9 Readers should be cautioned that the mass resolution of the reported 
simulation is 84000 M® = 5 M Je an S (125 cm~ 3 ,200 K). Ideally, if one prop- 
erly combines the refinement strategy and the molecular formation criteria, 
the local Jeans mass can be resolved this way without explicitly requiring it. 
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in our simulation represents a star-forming molecular cloud 
that is self-gravitating, is thus decoupled from the gas on 
the grid. 10 It is tagged with its mass Mmc, dynamical time 
fdyn =max([37r/ (32Gp gas )] 1//2 , 1.0 Myr), creation time f cr , and 
metallicity. Each molecular cloud particle gradually yields 
an actual stellar mass, M star (f), over 12 f<j yn which then con- 
tributes to stellar feedback (See Figure[T]and £|2.51 >. 

2.5. Stellar Feedback 

Observational evidence suggests that only ~ 2% of the 
gas in GMCs is converted into stars per dy namical time 
ilKrumholz & McKedl2005l: iKrumholz & Tanfeool and ref- 
erences therein). Numerical studies also indicate that turbu- 
lence, magnetic fields, or radiation pressu re can make the star 
formation process surprisingly slow (e.g. lMurray et al.ll2010l : 
I Wang et afll2010l) . To reflect these observations in our simu- 
lation, only 20% of the molecular cloud particle mass, Mmc> 
turns into an actual stellar mass, M stal (f ), over 12 ?d yn by 

M star (f)=0.2M M c / x'e- x 'dx' (6) 
Jo 

= 0.2M MC [1-(1+T) e - T ], (7) 

where T = (t — f C r)Adyn- In this formulation, the production 
of the stellar mass peaks at fd yn . As 7.5 x 10~ 7 of the rest 
mass energy of M star is gradually deposited into the cell in 
which the particle resides, 11 this thermal stellar feedback re- 
plenishes the energy loss to radiative cooling. At the same 
time, the rest of the molecular cloud particle mass, 0.8 Mmc, 
slowly returns to the gas grid. This again reflects the fact that 
most of the gas in GMCs does not end up locked in stars in 
a few dynamical time, but is blown out into the ISM to be 
recycled. Meanwhile, 2% of the ejected mass is counted as 
metals, contributing to the metal enrichment of the ISM (See 
Figure [TJ. 

Overall, our feedback treatment corresponds to the en- 
ergy of 10 51 ergs for every 750 M Q of actual stellar mass 
formed. Althou gh Type II supernovae explosions are it s dom- 
inant source (ISpitzerill990l:lTasker & Bry an 2006, 2008}, this 
feedback also models various other types such as pro tostel- 
lar outflows dLi & Nakamural l2006t iNakamura &~U[ 2008). 
photo ionization dMcKedl 19891) . and stellar winds dOev et al.l 
2001). Therefore no explicit time delay is necessary between 
the formation of a molecular cloud and the start of stellar feed- 
back. This thermal feedback heats the mass of - 10 4 M e in a 
< 30 pc cell up to ~ 10 K. but a multiphase medium_£Mc- 
Kee & Ostriker Tl977h is naturally established without using 
any su b-resolution model. The so-called overcooling prob- 
lem dSomerville & Primack|[T999l: iBalogh et al.l l2001) is ab- 
sent in our simulation since the cooling time of these hot cell s 
is much longer than the sound crossing time dKim etal.l2009l) . 

10 We point out that the usual terminology of star particle to represent 
10 5 - IU".W. has been a misnomer. We therefore make each of our particles 
to be 8000 Mq , regarding it as a molecular cloud gradually spawning stellar 
mass in it. These particles are still collisionless and do not fully represent the 
real nature of molecular clouds. However, we emphasize that our molecular 
cloud particles harbor a slow star formation rate matching observations. 

" Assuming the Salpeter initial mass function dn/dM °= (M/Mq) -2 3 in a 
star cluster (Salpeter 1955), the fractional mass which ends as Type II super- 
nova (SNII, > 9 M ) is 1.2%. Thus, fixing the mass of each SNII to be 9 M® 
we inject 10 51 ergs per 9 Mq/1.2% = 750 Mq of the stellar mass formed. This 
ratio lO 51 ergs/75OM = 1.3 x 10 48 ergsM ' equals to 7.5 x 10~ 7 of the 
stellar rest mass energy. 



2.6. Accreting Massive Black Hole 

A 1O 5 M massive black hole (MBH) is put as a seed at the 
center of each simulated galaxy. It is treated as a collision- 
less sink particle, but grows in mass by accreting gas from its 
surroundings. We estimate the rate of accretion by employing 
the Eddington-limited spherical Bondi-Hovle formula (Bondi 
& Hovle ll944tlBondilll952h : " 

M B h = min(M B , M Edd ) (8) 

. /47tG 2 M 2 H p B 47zGM BH m p \ 

= min — , — £ , (9) 

V cj e v a T c J 

where Mbh is the mass of a MBH, c s is the sound speed of the 
gas at the cell the MBH resides in, m p is the mass of a proton, 
and Or is the Thomson scattering cross-section. Note that, 
when compared with Eq.(Q3 the nondimensional parameter a 
is absent, pe is the density at the Bondi radius 

2GM BH ( M BH \ ^ 10km/s ^ 2 

RB= ^r c " 8 - 6pc {^mJ{^^) ,( 0) 

and is extrapolated from the density p ias of the cell of size Ax 
where the MBH resides by 

Pb =Pgas • min((A JC /7? B ) 1 - 5 ; 1.0) <p gas . (11) 

Here a n r~ 3 / 2 density pr ofile is assumed inside the sphere 
of R% dWang et alj l20ldh . Adopting a radiative efficiency 
£ r = 0.1 for a non-rotating Schwarzschild black hole (Shakura 
& Sunvaev TT973l:lBooth & Schavel 2009), the Eddington rate 
for a 10 5 M G black hole is ~ 0.002M e yr- 1 . For more dis- 
cussion on the accretion rate adopted here, see Appendix [Al 
To minimize any numerical artifacts, the gas mass accreting 
onto the MBH is uniformly subtracted from grid cells within 
a Bondi radius. The MBH also inherits the momentum of the 
accreting gas. 

Most importantly, to probe the gas dynamics accreting onto 
the MBH and to fully incorporate the MBH in a galactic sim- 
ulation, it is imperative to always reach the resolution close to 
the Bondi radius around the MBH. To resolve the gas around 
the MBH with the best resolution available, eight nearby cells 
close to the MBH are required to successively refine down 
to 15.2 pc (proper) at all times. In practice, the MBH natu- 
rally sits at the densest region most of the time, surrounded by 
many finest cells. While our spatial resolution is still slightly 
too large to resolve the Bondi radius of a 10 5 M Q black hole, 
Eq.dlOb. it is enough to resolve the Bondi radii of more mas- 
sive MBHs such as in nearby X-ray luminou s galaxies (e.g. 
- 120pc for SMBH in M87: lAllen eFaIll2006h . This shows 
that our simulations are beginning to depict the self-consistent 
coevolution of both galaxies and MBHs in one comprehen- 
sive framework. Admittedly, this resolution is still far from 
the Schwarzschild radius of any black hole 

2GM BH ,„_ 8 ( M BH \ „„ 

which is needed to thoroughly describe its accretion disk. Due 
to our resolution limit, a MBH particle in our framework rep- 
resents not just the black hole itself, but also includes accret- 
ing gas and stars deep within the galactic nucleus; in other 
words, the Bondi-Hoyle accretion estimate does not accu- 
rately model the physics below the resolution limit (See 35.21 . 
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| (A) radiative feedback | 
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| (B) mechanical feedback | 
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FIG . 2 . — Two dimensional schematic views of the different modes of massive black hole feedback. (A) radiativ e fee dback model described in £|2."7I photon rays 
carrying the energy are adaptively traced via full radiative transfer, (B) mechanical feedback model described in £|2.8| a momentum is injected to the cells around 
the MBH along pre-calculated directions, and (C) thermal feedback model predominantly used in particle-based galactic scale simulations: thermal energy is 
kernel-weighted to the neighboring gas particles around the MBH. 



2.7. MBH Radiative Feedback 

We now turn our attention to the feedback of an accreting 
massive black hole. The gravitational potential energy of the 
gas accreting onto a black hole is extracted during the gravita- 
tional infall. Assuming an infall down to the innermost stable 
orbit of an accretion disk, the conversion rate from the rest 
mass energy to feedback energy is 10%, previously defined 
as the radiative efficiency e r . Hence the bolometric radiation 
luminosity of a MBH is 

„2 



L BH = £rM m c 



(13) 



As was discussed earlier, for a long time a thermal energy de- 
position has been the dominant strategy to treat the feedback 
of an accreting MBH dSpringel et al.ll2005bt iDi Matteo et all 
l2005HColberg & Pi Matteol2008l:lBooth & Schavd20q9l: Cal- 
le garietal 12010 : Tevssie r et al.l2010L See lSijacki et al.l(l2008h 
or lDebuhr et alj (1201 ll) for other approaches). Without ques- 
tion, it has been an effective approximation characterizing the 
impact of an accreting MBH on a resolved scale when suffi- 
cient resolut ion or full radiative tra nsfer is inaccessible (See 
Figure 0C); Springel et al. 2005b). Despite its practical ef- 
ficiency, however, better feedback models are imperative for 
high-resolution galaxy formation studies where the Bondi ra- 
dius is starting to be resolved. In the next two sections, we 
explain the detailed implementations of two modes of MBH 
feedback: radiative and mechanical. The thermal feedback 
model previously used can be regarded as an approximation 
of these two feedback channels combined. 

Although the radiation from the MBH in a galaxy was 
tested in spherically symmetric or axisymmetric models 
(ICiotti & Ostriker 2007; Proga et a l. 2008; Ciotti et alf2009t 
iKurosawa et alJl2009b IPark & Ricottill2010h . a three dimen- 
sional radiative transfer calculation of the impact of a MBH 
has never been performed in galactic scale simulations. In 
what follows, we treat the MBH as a point source of radi- 
ation and carry out a three dimensional transport computa- 
tion to evolve the radiation fields (See Figure |2j A); note that 
molecular cloud particles are not treated as radiation sources). 
Achieving high resolution around the MBH is critical here be- 
cause, if otherwise, the optical depth of the radiation could be 
small even at the s mallest resolved distance from the MBH 
(lOmma et al.ll2004l) . 

Enzo's radiative tra nsfer module incorpora tes the adaptive 
ray tracing technique (lAbel & Wand elt 2002) with the hydro- 



dynamics, energy, and chemistry solvers. It has been applied 
to problems such as the radiative feedback from Pop III stars 
(lAbel et al.ll2007t IWise & Abel 120081: 1 Wise et al.ll2010h and 
from Pop III black holes dAlvarez et al. I I2009D . For algorith- 
mic and numeri cal details of Enzo radiative transfer we refer 
the readers to IWise & Abell (H010); and, here we briefly de- 
scribe the machinery relevant to the presented results. First 
the luminosity of the MBH is assigned by Eq.(fT3l). Then 768 
(— 12 x 4 3 ; Healpix level 3) rays are isotropically cast with a 
monochromatic energy of £p n = 2 keV, a characteristic tem- 
perature of an averaged quasar spectral energy distribution 
(SED: ISazonov et al.ll200l. 120051 : ICiotti & Ostriker|[2007h . 12 
Consequently the number of photons per each initial ray is 



^init — 



dt ph _ £ r (M BH <ifph)c 



-'ph 



768 



£ ph -768 



(14) 



given the photon timestep dt p h which we set as the the light- 
crossing time of the entire computational domain. This choice 
is justified because the photons are in a free streaming regime, 
and the energy deposited by the radiation per timestep is rela- 
tively small. Each ray is traced at speed c until the ray reaches 
the edge of the computational domain or most of its photons 
(99.99995%) are absorbed. It is adaptively split into four child 
rays whenever the area associated with a ray becomes larger 
than 0.2 (Ax) 2 of a local cell. 

Photons in the emitted ray then interact with the surround- 
ing gas in three ways: they (1) ionize the gas, (2) heat the 
gas, and (3) exert momentum onto the gas. First, the ray loses 
its photons when it photoionizes H, He, and He + with the re- 
spective photoionization rates of 

Pjl-e-^)(E ph Y ktH /E iM ) 



£ph,H = 



^ph.He : 



n H (Ax) 3 dt ph 
-e-^)(E ph Y k . He /E, Me ) 



n He (Ax) 3 dt ph 



Mi 



T He+ N 



n He + (Ax) 3 dt p h 



(15) 



(16) 



(17) 



12 A characteristic temperature of the quasar SED is estim ated by equat- 
ing Compton heating and Compton cooling by the given SED I Sazo nov et all 
120041) . Therefore it can be considered as the temperature of a Comptonized 
hot plasma in the vicinity of the MBH, which is represented by our MBH 
particle resolved only by 15.2 pc. Note, however, that the choice of the 
monochromatic photon energy, £ p h, does not change the total luminosity of 
the MBH, nor does it greatly affect the nonrelativistic Klein-Nishina cross 
section for Compton scattering, o"kn, m me regime of £ p h <C m e ir. 
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TABLE 1 
Simulation Suite Description 



Physics a 


Sim-SF 


Sim-RF 


Sim-MF 


Sim-RMF 


Molecular cloud formation 


(See %2A\ 


o 


o 


o 


o 


Stellar feedback 


(See 323) 


o 


o 


o 


o 


Massive black hole accretion 


(See [|2l>) 


o 


o 


o 


o 


Massive black hole radiative feedback 


(See ED 


X 


o 


X 


o 


Massive black hole mechanical feedback 


(See 


X 


X 


o 


o 



n Foi~ detailed explanation, see the referenced section, o = included, x = not included. 



where P m is the number of photons coming into the cell, 
= nftGftdl is the optical depth, «h is the hydrogen num- 
ber density, Oh is th e energy-dependen t hydrogen photoion- 
ization cross-section ( iVerner et alJI 19961) , dl is the path length 
through the cell, and Ej = 13.6, 24.6, 54.4eV are the ioniza- 
tion thresholds for H, He, He + , respectively. The factors Y^ 
are the energy fractions used for ionization when secondary 
ionizations are considered (Shull & van Steenberdl985l) . 13 

Second, the excess energy above the ionization threshold, 
Ei, heats each of the species with the photoheating rates of 

r H = - 7T ' P - , etc. (18) 
nu(Ax) 3 dt vh 

where is the fraction of energy deposited as heat when sec- 
ondary ionizations are taken into account. 14 The 2 keV soft 
X-ray photon can also scatter off and heat an electron result- 
ing in the Compton heating rate of 

_ p in (i- e -*)AE(r e ) 

ic ~ 71 \TT. (-19) 

n e (Axydt ph 

where T e = n c Ga^dl is the optical depth, « e is the elec- 
tron number density, Okn is the nonrelativistic K lein- 
Nishina cross-section (~ Or ; Rybicki & Lightman il979h . and 
AE(T e ) = 4k B T e ■ (£ph /m e c 2 ) is t he nonrelativistically tra ns- 
ferred energy to an electron at T e dCiotti & O striker 200]]) • It 
should be noted that, in Compton scattering, a photon loses 
its energy by a factor of AE (T e ) / E^, but essentially keeps 
propagating without being absorbed. However, in order to 
model this with monochromatic photons, we instead subtract 
Pj n (l — e _Tc )A£ , (r e )/£ , ph photons from the ray. This is an- 
other way a ray loses its photons while traveling through a 
cell. Combined, the total heating rate by absorbed and scat- 
tered photons becomes 

T = "Hr H + «H e r H e + "He+ r He+ + n e r C • (20) 

Lastly, photons exert outward momentum to the gas when 
they are taken out from the ray either by photoionization or by 
Compton scattering. It was claimed that the radiation pres- 
sure from the MBH may markedly alter the environ ment near 
the MBH, especially withi n ~ 0.1 kpc in radius (lHaehneltl 
Il995t iDebuhr et al.l 1201 ll) . The large-scale galactic wind 
driven by deposited photon momentum is also considered as 
a possible explanation for the Mrh — Oh,.i 2f . relation (Murray 
et al. l2005l) . The added acceleration onto the cell by the radi- 
ation pressure is calculated by 

fl'Pph Aost^ph j. 

ph m ce \idt p h p g!is {Ax) 3 cdt ph 



(21) 



13 Y kM =0.3908(1 -#4092)1.7592 and y tHe = .0554(1 -* 04614 ) L66S0 are 
fitted as a function of an ionization fraction x = n H + /fiH,iot — «He+ / n He tot I 
the effect of secondary ionizations on He + can be ignored. 

14 Note that Y T = 0.9971 [1 - (1 2263)1.3163] approac hes when the 
ionization fraction gets close to 0. In other words, when the ionization frac- 
tion is low photons are preferentially used to first ionize the gas rather than to 
heat the gas. 



where t/p p h is the photon momentum exerted onto the cell in 
dt p h, ^lost is the number of photons lost in the cell, and f is 
the directional unit vector of the ray. Neglecting the radiation 
pressure on dust grains is conservative because its inclusion 
would further enhance the negative feedback effect (See £]5.2| ). 

2.8. MBH Mechanical Feedback 

Observations find that a significant portion of the energy 
extracted during the accretion onto a M BH is released as me 



chanical energy, cre ating bipolar jets (Bridle & Perlevl 1984 : 
Pounds et al. 2003) or inflating cavities dFabian et al.ll2002t : 
iMcNamara et al.ll2005l) at the sites of active galactic nuclei 
(AGN). A number of authors have used a numerical approach 
to explore the effect i veness of jets in heating u p a cooling flow 
dFabian et alJll994l ; |Peterson & Fabianll2006h ; most of them 
targeted the gas dynamics in galaxy clusters with ~ kpc res- 
olution excluding detailed galactic scale phvsics (e.g. Omma 
et al. 2004; Cattaneo & Tevssier 2007; Antonuccio-Delogu & 
Silk l2008t Dub ois et alj|2010l) . In the meantime, a numerical 
analysis on stellar winds from nuclear disk or MBH jets has 
been carried ou t in a galactic scale, but only in an one d imen- 
sional context dCiotti et al.l2009ll2010fc[Shin et alj2010|) Here, 
we construct a mechanical feedback model of a MBH appli- 
cable in three dimensional galactic simulations, which creates 
accretion-rate-dependent subrelativistic bipolar jets launched 
at the vicinity of the MBH (See Figure^B)). 

Let us assume that all of the bolometric luminosity of the 
MBH, Lbh, is converted to the "mechanical" power of jets. 
Because the ejecta has to climb out of the potential well of the 
MBH, the "kinetic" power of the jets is less than Lbh by 

^BH = ^mech (22) 

= fldn + (gravitational potential energy). (23) 

Therefore the "kinetic" power of the jets, as we introduce at a 
scale of 7?j e t = 2Ajc = 30.4 pc, can be written as 

fkin = ekin^BH = S^E.MnnC 1 = -MjetV? t , (24) 

where g^in < 1 is the "kinetic" coupling constant denoting the 
fractional energy available for the kinetic motion of the jets 
(See Figure |2fB)). Mj e t is the mass ejection rate of the jets, 
and vjet is the jets velocity when introduced in the simulation. 
Hence £ki n encapsulates not only the acceleration of the jets 
powered by the AGN central engine, but also the gravitational 
"redshift" from the scale of an accretion di sk (~ R<j c h) t o a re- 
solved scale of jets in simulations (~ Rj* ). [Ciotti et all (120091) 
provides estimates for a MBH of / = M BH /M E dd = 0.005 as 

s 0.0015 (25) 



Tict 



e T M BH c 2 e r (l+4000 4 
0.2 



_ M jet 



M BH (1 + 100/) 



0.04, 



(26) 
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FIG. 3. — A projected density of the simulation box (16 comoving Mpc) at z = 3 is displayed on the right; circles represent the identified massive halos. On the 
left a 9.2 x 1O"M0 halo, i.e. the model galaxy, is shown in a 200 kpc box (proper). High-resolution images are at http://www.jihoonkim.org/. 



based on which we fiducially adopt conservative values of 
£km = 10~ 4 and r/jet = 0.05. 

With e^in and r/j et now fixed, the kinetic motion of the jets 
can be fully described. First, as usual, out of a sphere of Rg 
centered on the MBH the accreting mass is taken out at every 
finest hydrodynamical timestep dt; then 5% of the accreted 
mass, M jel dt = 0.05 M^ndt, is set aside as a mass of jets. Now 
Eq. (l24b yields the initial jet momentum, (Mj et c/f )vj e t, with 



Vjet 



V fJjet 



1/2 



= 6000 km s~ 



(27) 



for £ r = 0.1. This value o f Vj e t is consistent with numerous 
observational evidence (e. g . iBiretta & Junorl 19951 : lJunor et ail 
119991: iHoman et alJl2009h and relat ivistic MHD simulations 
(e.g. lVlahakis & K5n igl 2003, 2004) suggesting the existence 
of at least mildly relativistic AGN jets on scales of 1 - 10 pc 
from a central engine. This Vj et is also well-matched wi th the 
veloci ty of momentum-driven AGN winds discussed by lKingl 
(2009). Finally the launch speed of the surrounding cells is 
found by averaging t he momentum of jets and the preexisting 
gas in those cells (IWang et al.ll2010h . 

One may want to continuously launch the jets at every finest 
hydrodynamical timestep. However, if the injected mass of 
jets is minuscule compared to the preexisting mass in sur- 
rounding cells, the jets make little or no dynamical impact 
on the surrounding cells after being mass-weighted averaged 
with them. Since it is unfeasible to resolve all gas cells around 
the MBH down to Mj ct dt, an alternative approach is indispens- 
able. Moreover, there is growing observational evidence of 
double-lobed radio galaxies (or double-double radio galaxies; 
DDRG) implying that the jets have launched in an episodic 
fashion with jets interruption tim escales of 10 5 -10 8 years 
(IStawarzll2004l : ISaikia et al.ll2006h . These two considerations 
lead us to adopt the following method: every time the accu- 
mulated jet mass, T.Mj et dt, exceeds the threshold of 300 M 
it is injected in collimated bipolar jets of a width of five finest 



cells in the vicinity of the MBH. This approach renders jets 
intermittent (once every 30 Myr if Mbh = 10 _5 M Q yr _1 ) and 
dynamically important in our calculation. 

The jets are injected parallel and anti-parallel to the total 
angular momentum L of the accreted gas up to that point. The 
angular momentum vector L changes its direction frequently 
while it asymptotes to the overall galactic rotation axis. This 
implementation is motivated by the observations of X-shaped 
radio galaxies (XRGs) where the radio jets rapidly reorient 
themselv es by the interaction w i th the surrounding gas or b y 
mergers (Merritt & Ekers1 l2002l: fGopal-Krishna et al.ll2003h . 
Lastly, since the mechanical or the radiative feedback alone 
may not describe the whole picture, we include hybrid models 
in which each of these two channels constitutes half of the 
MBH bolometric luminosity, Lbh (Sim-RMF; see Table[T|i. 

Note that the mechanical channel has not been a main driver 
of MBH feedback in the presented calculation because, with 
highly suppressed mass accretion rate, jets have launched only 
a few tens of times in 350 Myr (See §4.2\ , We later comment 
upon its efficiency in 



3. INITIAL CONDITIONS 

These improved physics of galaxy formation are first ex- 
tensively tested in isolated galaxies. We then apply them to a 
massive star-forming galactic halo of 9.2 x 10 Mq at redshift 
3 in a cosmological ACDM simulation. We begin by describ- 
ing how the initial conditions of our simulation are generated. 



3.1. 



Setting up a ~ 10 Mq halo 

A three dimensional cubic volume of 16 comoving Mpc 
on a side is set up at z = 60 assuming a flat ACDM cos- 
mology with dark energy density Q.a = 0.73, matter density 
i2 m = 0.27, baryon density ilj, = 0.044, and Hubble constant 
h = 0.71 (in the unit of H = 100km s _1 Mpc _1 ). A scale- 
invariant primordial power spectrum (spectral index n = 1, 
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lEisenstein & Hdfl999l) is adopted with <7 8 = 0.81, the rms 
density fluctuation amplitude in the sphere of 8 /z^'Mpc. 

We identify a dark matter halo of ~ 1O 12 M at z = 3 by per- 
forming a coarse-resolution adiabatic run. Then we recenter 
the density field around this halo and set up a new initial con- 
dition which preserves the same large-scale power yet con- 
tains a small-scale power as well, with a 128 3 root grid and a 
series of two nested child grids of twice finer resolution each 
(160 3 cells for level / = 1, and 200 3 for I = 2). Therefore the 
finest nested grid at level / = 2 spans 6.25 comoving Mpc on a 
side, contains 200 3 dark matter particles of 9.6 x 1O 5 M , and 
manifests the equivalent resolution of a 512 3 unigrid. Initially 
all the cells throughout I — 2 grids are allowed to be further re- 
fined; however, the volume in which additional refinement is 
enabled (V re f ; in the shape of a rectangular solid) continually 
shrinks in size in such a way that it encloses only the smallest 
dark matter particles. 15 An initial metallicity of Z = 0.003 Z Q 
is also set up everywhere to track the metallicity evolution and 
to facilitate cooling below 10 4 K. 

Our initial condition is first evolved to z = 3 with a low- 
resolution (121.6 pc) refinement strategy and a particle for- 
mation and feedback recipe, without an accreting MBH. At 
z = 3 we split each dark matter and star particle inside the 
focused volume (Vf oc ; a rectangular hexahedron of 1.28 co- 
moving Mpc on a side, a subset of V re f) into 13 child parti- 
cles using the particle refinement technique by Kitsionas & 
Whitworth (2002). This algorithm places child particles on a 
hexagonal close packed (HCP) array, and has been applied to 
many particle-based applications requ iring enhanced particle 
resolution in a resimulated region (e .g. Bromm & Loeb1 l2003l : 
iKitsionas & Whitworthll2007t I Yoshida et al.ll2008h . After the 
particle splitting procedure, each dark matter particle in Vf oc 
represents a collective mass of 74000 M . Across Vf oc cells 
are now allowed to refine up to 11 additional levels, achieving 
maximum spatial resolution of 15.2 pc at z ~ 3 (See £|2.21 i. 



3.2. Galactic parameters 

Consequently, this process produces our focused object at 
z = 3 dubbed a model galaxy, on which a suite of high- 
resolution simulations is performed (Figure |5J. The model 
galaxy has a mass of M v ; r ~ 9.2 x lO n M at z — 3 and a cor- 
responding virial radius of 



R v 



M, 



1/3 



2G 



-1/3 



310 comoving kpc(28) 



given = 0.71, Q. m = 0.27, and A c = 200. The dark mat- 
ter halo represented by ~ 1.1 x 10 7 particles constitutes ~ 
88% of the total mass. About ~ l.Ox 10 7 particles contain 
8.0 x 1O 1O M of stellar mass, whereas the rest, 3.5 x 1O 1O M , 
is in gaseous form available for future star formation, either 
in the ISM or in the embedding halo. There is no shortage 
of gas supply, as the gas from outside the halo continuously 
falls inward either by spherical accretion or by cold accre- 
tion along one of th e multiple filaments dDekel et alj|2009l : 
ICeverino et alj |2010). The halo has spin parameters of 0.051 
for dark matter, and 0.069 for gas. At the center of all lies a 
1O 5 M MBH we plant as a gravitati onal seed. This choice of 
the initial MBH mass lies below the Magor rian et alJ (Il998h 

15 This active adjustment on the size of V re f prevents heavier dark matter 
particles of initial 1 = and 1=1 grids from penetrating the central region of a 
simulation box, thereby causing runaway refinement. Typically V rc f becomes 
~ 60% of the entire 1 = 2 region in length at z = 3, which is still large enough 
to encompass the Lagrangian volume of a ~ 10 halo at z = 3. 



Sim-SF 
Best fit 
KennicutK1998) 




log 2 aas (M SU nP c ' ) 



-gas 

FIG. 4. — The relationship between star formation rate (SFR) and gas sur- 
face density. The data is from a (20 kpc) 3 box centered on a MBH in Sim-SF 
at z = 2.75. The solid line is the best fit for simulated dat a while the dash ed 
line is from observations of nearby galaxies (Esfr x ^las^ Kennicutt 1998). 

relationship assuming 10% of the stellar mass is in the bulge, 
which may have resulted in a weaker mode of MBH feedback 
- possibly a "radio-mode" a nalog ue - and the negligible gas 
expulsion by the MBH (See $Q1 and Sj4~4l , Therefore the re- 
ported results should not be interpreted as a general picture of 
MBH feedback. It remains to be seen whether more massive 
MBHs or fast growing MB Hs h ave different effects. We will 
come back to this issue in ^4.41 



4. RESULTS 

A suite of simulations with optional modes of feedback is 
performed from z = 3 to 2.6 in order to investigate the evolu- 
tion of a massive star-forming galaxy with its embedded mas- 
sive black hole. We mostly focus on two simulations, one with 
and the other without MBH feedback (Sim-SF and Sim-RMF; 
TableHJ. Each of the calculations is performed on 16 proces- 
sors of the Orange cluster 16 at Stanford University. Grids and 
particles altogether, each simulation is routinely resolved with 
~ 6.5 x 10 7 total computational elements (~ 4.5 x 10 7 parti- 
cles and ~ 270 3 cells). To evolve the system for 350 Myr, 
each of these runs typically takes ~ 20000 CPU hours. 

4. 1 . Star Formation Rates 

First we ch eck t he validity of our mol ecula r cloud forma- 
tion criteria ( ^2.4b and stellar feedback (i ]2.5l ) by comparing 
star formation rate (SFR) with gas density. Figure [4] displays 
a relation between the SFR surface density and the gas surface 
density in Sim-SF at z = 2.75. In a (20 kpc) 3 box centered on 
a MBH, each data point is made by taking the mean values in 
a (1 kpc) 3 bin, the typical aperture size of £ spr — E gas stud- 
ies fo r spatially-resolved nearby galaxies (e.g. Kennicutt et al. 
l2007h . Here Eq.(|7} is used to calculate the stellar mass newly 
spawning in each cell, and the data points below the observa- 
tion limit, lO _5 M yr _1 kpc -2 , are discarded. 

The molecular cloud formation and stellar feedback, joined 
with high spatial resolution, work together to self -regulate star 
formation. However, authors note that our best fit to this par- 
ticular snapshot of the galaxy is steeper than the observed 
trend of z ~ with larger dispersion. 

16 Infiniband-connected AMD, 8 cores per node, 4 GB memory per core 
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FIG. 5. — The face-on views of the disks. Density in the central 20 kpc (proper) sliced through the MBH (black dots at the centers) at z = 2.75, about 220 Myr 
after the MBH is placed. Sim-SF on the left, and Sim-RMF on the right. Compare with Figure[6] 
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FIG. 6. — The face-on views of the disks. Temperature in the central 20 kpc (proper) sliced through the MBH (white dots at the centers) at z = 2.75. Sim-SF 
on the left, and Sim-RMF on the right. A hot region of a size ~ 2 kpc in Sim-RMF heated by MBH feedback is prominent, which remarkably contrasts with a 
much colder ISM in Sim-SF. 
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FIG. 7. — The mass-weighted radial profile of temperature in a 20 kpc 
sphere centered on the MBH at z. = 2.75. The red solid line and the green 
dashed line represent Sim-SF and Sim-RMF, respectively. The temperature 
within ~ 2 kpc radius is raised mostly by the radiation from the MBH. 



4.2. Lack of Star- forming Gas in the Inner Core 

Now we turn to the topic of an accreting massive black hole 
and its feedback. We focus on how MBH feedback changes 
its surrounding ISM, and how it locally suppresses molecular 
cloud formation. For this purpose, we hereafter examine the 
snapshot of the model galaxy at z = 2.75 (or 2410 Myr af- 
ter the Big Bang), about 220 Myr after an accreting MBH is 
placed at the center of the galaxy. The model galaxy now has 
a mass of M v ; r ~ 9.0 x 10 M and correspondingly, a virial 
radius of R v j r ~ 80 kpc (proper; radii are hereafter in proper 
kpc, not comoving, unless marked otherwise). 

When a MBH starts to accrete gas, the gravitational po- 
tential energy of the accreting gas is released in the form 
of radiation and jets. Even in the case of a slowly growing 
MBH, as in our simulations (a possible "radio-mode" ana- 
logue; M B h ~ 0.05M Edd ; see £|4.4| ), the feedback from the 
MBH is known to play a major role in regulat ing star forma- 
tion and its own growth (Croton et al. 2006; McNamara & 
Nulsen l2007tlSiiacki et al.ll2007l ; lPuchwein et al.ll2008h . 
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FIG. 8. — Joint probability distribution functions (PDFs) of gas density and temperature colored by gas mass in each bin. The data is for a 200 kpc sphere 
centered on the MBH at z = 2.75. Sim-SF on the left, and Sim-RMF on the right. The vertical dashed line in each plot denotes the density threshold for molecular 
cloud formation « t hres ( fl2.4t . Note that the MBH feedback in Sim-RMF heats the dense gas up to 10 6 -10 7 K increasing the amount of high density gas stable 
against fragmentation (zone "A").. 



The density and temperature structures in the central re- 
gions of the galaxies from Sim-SF (left; without MBH feed- 
back) and Sim-RMF (right; with radiative and mechanical 
MBH feedback) are shown in Figures [5] and [6] In particular, 
in Figure [6] for Sim-RMF, a hot region of size ~ 2 kpc sur- 
rounding the MBH is prominent, which remarkably contrasts 
with a much colder ISM in Sim-SF. This region is heated up 
to 10 6 K by ionizing photons heating hydrogen and helium, 
and scattering off electrons. The latter, i.e. Compton heating, 
is important especially in a highly ionized region. When the 
fractions of neutral species (or singly ionized He) are low, the 
effect of photoheating on H, He, and He + is mild; instead, 
the contribution of Compton heating on electrons is relatively 
large. The hot temperature at the center of Sim-RMF is also 
evident in the radially averaged temperature profile of Figure 
[7] Note that the 1-2 kpc distance over which the gas is heated 
is consistent with the c haracteristic distance found in an ana- 
lytic study (Figure 4 of ISazonov"eT al. 2005) out to which the 
gas is heated by photoionization and Compton scattering. 

A hot temperature in the inner core of the galaxy in Sim- 
RMF leads to a significant deprivation of cold, dense star- 
forming gas. Figure [8] illustrates how the structure of ISM 
is changed by MBH feedback, in terms of joint probability 
distribution functions of gas density and temperature. 

• The left figure depicts a typical ISM without MBH 
feedback but still with stellar feedback (Sim-SF). It 
features a multiphase ISM that is naturally achieved 
in adaptively refined mesh, including cold, dense star- 
forming gas (T < 10 4 K, p > 10~ 24 gcm~ 3 ), and 
hot diffuse supernovae bubbles. As expected, above 
the molecular cloud formation threshold (n tn i- e s = 125 
cm~ 3 ; denoted by a dashed line), gas cells immediately 
turn into molecular cloud particles, and thus no cell is 
left behind unresolved. 

• On the right, the gas cells of density > 10~ 24 g cm~ 3 
are now heated up to 10 6 - 10 7 K, populating zone "A". 



These cells are mostly located on the disk in relatively 
close proximity (< 2 kpc) to the central MBH. These 
cells are very stable against fragmentation because they 
are so hot that they can barely cool down to a typical 
molecular cloud temperature in a dynamical time (i.e. 
fcool S> fdyn)- For that reason, the cells have hard time 
to fulfill the condition (c) of the molecular cloud forma- 
tion criteria described in ^2A\ The heating by the X-ray 
radiation thus increases the amount of dense gas in the 
vicinity of the MBH which is incapable of condensing 
into stars. 17 

Therefore, the feedback from even a slowly growing MBH 
retains hot dense gas at a galactic center which otherwise 
could have created strong star formation. Figure [9] dramati- 
cally demonstrates the distinct changes in radially averaged 
density profiles when MBH feedback is included. 

• The left figure of gas density profiles displays that the 
radiation from the MBH keeps a substantial amount of 
gas at the core of the galaxy, and inhibits the gas from 
all turning into stars. The total gas density of Sim-RMF 
(the green dashed line) at 0.1 kpc from the MBH is 
about ~ 5 times as high as that of Sim-SF (the red solid 
line). The thin blue dotted line represents the initial pro- 
file of gas at z=3 when the simulation restarts with 15 
pc resolution. However in Sim-RMF, only a minimal 
fraction of the gas within 0.2 kpc radius is below 10 4 K 
(the cyan dot-dashed line) whereas in Sim-SF, 10 - 30% 
of the gas in the same region is considered to be cold, 
thus potentially star-forming (the pink dotted line). 

• The deprivation of cold dense gas in Sim-RMF in- 

17 Some hot gas cells above the molecular cloud formation threshold still 
exist on this PDF, not turning into particles (zone "B"). While it is due to the 
violation of one of the molecular cloud formation conditions (f coo i < fd yn ), we 
emphasize that the thermodynamical properties here are unreliable as they are 
above the resolution-dependent molecular cloud formation threshold. 



Galaxy Formation with Stars and MBHs 



11 



10" 



Cl 



10° 



10' 



o 

§9 s 



10* 



Sim-SF-gas 
Sim-RMF-gas 
baseline at z=3 
Sim-SF-gas-cold 
Sim-RMF-qas-cold 



■'■"■"'•X. 



10" 



Sim-SF 
iim-RMF 



0.1 



Radius (kpc) 



10 




Radius (kpc) 




Sim-SF-gas 
Sim-RMF-gas 
Sim-SF-gas-cold 
Sim-RMF-qas-cold 



FIG. 9. — Spherically-averaged radial gas density profiles centered on the MBH at z = 2.75. In each panel, the red solid line and the green dashed line represent 
Sim-SF and Sim-RMF, respectively. Left: in Sim-RMF the radiation from the MBH keeps a large amount of gas at the galactic center (green dashed line) 
compared to Sim-SF (red solid line). Only a minimal fraction of the gas within 0.2 kpc radius is below 10 4 K (the cyan dot-dashed line) whereas in Sim-SF, 10 - 
30% of the gas in the same region is considered to be cold (the pink dotted line). The blue thin dotted line shows the density profile at z = 3. Right: star formation 
rate (SFR) density shows significantly suppressed star formation activity in the central 0.2 kpc sphere. 

within 0.2 kpc than what Sim-SF does (5.1 x 1O 7 M ), but al- 
most none of this gas is below 10 4 K. This gas in the core is 
not consumed by star formation, nor is pushed away by any 
mechanical outflow. Note that in Sim-RMF cold gas mass 
within 10 kpc is reduced by ~ 50% displaying how far the 
MBH radiation reaches. 

Readers should also note that the enclosed gas masses at 
the virial radius (80 kpc proper) are almost identical between 
Sim-SF and Sim-RMF, implying there has been no massive 
gas expulsion driven by the MBH. This is one of the key differ- 
ences from previous numerical studies. Very strong gas expul- 
sions were frequently observed in previous studies in which 
the MBH feedback energy is only thermally deposited to a 
few neighboring gas particles (e.g. Spr ingel et al.ll2005bl) . In 
contrast, most of the gas is still bound to our simulated galax- 
ies because (1) the energy released from our slowly grow- 
ing MBH is relatively small (a possible "radio-mode" ana- 
logue), 19 (2) the gas mass to which the MBH energy is cou- 
pled is large in our radiative feedback formalism, and (3) the 
mass accretion rate onto our MBH is not high enough to re- 
peatedly drive jets. Note again that the mechanical channel 
of MBH feedback is not a main driver of feedback in the pre- 
sented calculation. The MBH has not double d its mass at the 
end of our calculation (after 350 Myrs; see ^4.4b ; and, with 
this highly suppressed mass accretion rate, jets have launched 
only a few tens of times in 350 Myr. We com e back to the 
efficiency issue of mechanical feedback in £15.21 

To summarize, we have shown that MBH feedback, espe- 
cially its radiation, alters the multiphase ISM of the surround- 
ing gas and thus deprives the galactic inner core of cold, dense 
star-forming gas. Two consequences arise from the lack of 
star-forming gas at the galactic center: locally suppressed star 
formation, and the associated change in stellar distribution. 
We discuss these topics in the following sections. 

9 This is valid only for the results presented herein in which the mass ac- 
cretion onto the MBH is highly suppressed by self-regulation. It remains yet 
to be seen whether more massive MBHs or fast growing MBHs (for example, 
in merging galaxies) do or do not unbind the gas from the galactic gravita- 
tional potential. T he so -called "quasar-mode" feedback will be the topic of a 
future paper (See $5.2\ . 
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FIG. 10. — Enclosed mass profiles of total gas and cold gas (T < 10 4 K) 
centered on the MBH at z = 2.75. The red solid line and the green dashed 
line represent Sim-SF and Sim-RMF, respectively. Sim-RMF harbors more 
gas within 0.2 kpc than what Sim-SF does, but almost none of this gas is 
below 10 4 K. 

evitably prompts the suppression of star formation ac- 
tivity in the inner core of the galaxy. The right figure re- 
veals the star formation rate density (in M yr kpc" 3 ) 
as a function of distance from the MBH. Eq.(|7]l is again 
used to calculate the new stellar mass being generated 
in each cell. The SFR density of Sim-RMF at 0.1 kpc 
from the MBH is reduced by more than ~ 50% when 
compared with that of Sim- SF. 18 Overall, the X -ray ra- 
diation from the MBH severely suppresses the SFR of 
Sim-RMF inside the 0.2 kpc sphere. 

The gas masses enclosed within a radius r, M gas (< r), are 
plotted in Figure[l0] showing the impact of MBH feedback as 
a powerful energy source to reshape the galactic gas distribu- 
tion. Sim-RMF harbors ~ 2.5 times more gas (1.3 x 1O 8 M ) 

18 Note that since we used Eq.0 to estimate SFR, each molecular cloud 
particle will generate stellar mass for 12 fd yn . Therefore the SFR inside the 
0. 1 kpc sphere would include a number of particles that were formed outside 
the sphere, but have migrated inward and are now forming stars there. This is 
why the SFR density is not as suppressed as one would have naively expected 
from the density profile of cold gas. 
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FIG. 11. — Stellar density profiles of total stars and young stars (molecular 
cloud particles of age < 100 Myr) in a 20 kpc sphere centered on the MBH 
at z = 2.75. The red solid line and the green dashed line represent Sim-SF 
and Sim-RMF, respectively. Locally suppressed star formation at the galactic 
center in Sim-RMF leads to a considerable reduction of stellar density in the 
region. 
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FIG. 12. — Enclosed mass profiles of total stars and young stars (molecular 
cloud particles of age < 100 Myr) centered on the MBH at z — 2.75. The 
red solid line and the green dashed line represent Sim-SF and Sim-RMF, 
respectively. The stellar mass enclosed in the 0. 1 kpc sphere of Sim-RMF is 
about an order of magnitude smaller than that of Sim-SF. 

4.3. Locally Suppressed Star Formation and the Change in 
Stellar Distribution 

The inner core of the galaxy in Sim-RMF becomes a sterile 
environment for molecular cloud formation (star formation) 
because the gas is hot and turbulent, therefore Toomre stable. 
As a consequence, star formation is suppressed locally in the 
inner core, as shown in Figure|9] Figure[TT]displays how the 
stellar mass density profile changes in the inner core region as 
a result of the locally suppressed star formation in Sim-RMF. 
Again we use the snapshot at z — 2.75, about 220 Myr after 
the MBH is placed at the center of the model galaxy. The 
stellar mass density at 0. 1 kpc in Sim-RMF is only less than 
~ 20% of that of Sim-SF. The mass density of young stars 
(age < 100 Myr) shows the similar drastic reduction. 

The stellar masses enclosed within a radius r, M stal (< r), 
are shown in Figure [12] The stellar mass enclosed within the 
0. 1 kpc sphere of Sim-RMF is almost an order of magnitude 
smaller than that of Sim-SF. The suppressed star formation 
replaces the steep inner cusp of stellar density profile with 
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FIG. 13. — Star formation history in a 0.2 kpc sphere: the mass of new 
stars (molecular cloud particles born after the MBH is placed at 2190 Myr) 
inside a 0.2 kpc sphere centered on the MBH. The red solid line and the green 
dashed line represent Sim-SF and Sim-RMF, respectively. Since 2300 Myr 
the star formation activity of Sim-RMF in this region is suppressed. 

a flattened core ~ 0.3 kpc in radius; i.e. stars are less con- 
centrated at the galactic center of Sim-RMF. Considering the 
relatively small amount of energy released from the slowly 
growing MBH, the difference between these two lines is quite 
remarkable. Together with Figure [10] one expects that the 
stellar to gas mass ratio inside the 0.2 kpc sphere of Sim- 
RMF will be much smaller than that of Sim-SF. Note also 
that the total stellar mass enclosed at the virial radius (80 kpc) 
is almost indistinguishable between Sim-SF and Sim-RMF. 
In other words, star formation in Sim-RMF is not globally 
suppressed, but only locally suppressed at the center. This is 
because our MBH feedback is not strong enough to unbind a 
large amount of gas (See ^4.2t , or to globally abolish cold, 
star-forming clumps in the entire ISM. 

Figure Qj] exhibits the evolution of the new stellar mass 
(molecular cloud particles born after the MBH is placed at 
2190 Myr, or at z = 3) inside a 0.2 kpc sphere centered on the 
MBH. The plot demonstrates that since 2300 Myr the star for- 
mation activity in the inner core of Sim-RMF is suppressed. 
Naturally, alteration in the stellar distribution ensues at the 
center of the galaxy with an active MBH. Figure [14] strikingly 
contrasts the distribution of newly-formed stars (molecular 
cloud particles of age < 10 Myr). 

• In the top row, a three dimensional rendering of newly- 
formed particles is constructed at a ~ 45° angle from 
the disk plane. Here the light intensities of newly- 
formed particles are integrated along the lines of sight. 
At the center of the stellar distribution, i.e. the densest 
peak, lies the MBH particle. 

• In the bottom row, newly-formed stellar masses are pro- 
jected along the z-axis of the simulation box which 
makes a ~ 47.2° angle with the angular momentum 
vector of the gas within a 5 kpc sphere centered on the 
MBH. The morphological difference at the inner core of 
the galaxy is particularly evident. Stars in Sim-SF are 
highly concentrated at the center, while stars in Sim- 
RMF are less concentrated but form spiral-like struc- 
tures at ~ 0.2 kpc radius from the MBH. 

In summary, it is shown that MBH feedback suppresses 
star formation locally at the galactic inner core, thus signif- 
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FIG. 14. — The distribution of newly-formed stars (molecular cloud particles of age < 10 Myr) at z = 2.75. Sim-SF on the left, and Sim-RMF on the right. 
Top: images of newly-formed stars constructed at a ~ 45° angle from the disk plane. Visualization courtesy of Ralf Kaehler. Bottom: newly-formed stellar mass 
projected along the z-axis of the simulation box which makes a ~ 47.2° angle with the gas angular momentum vector. Each frame is 3.0 kpc wide. The difference 
in stellar distribution is dramatic, especially in the < 0.5 kpc core. 



icantly changing the stellar distribution there. This new chan- 
nel of feedback is particularly interesting because it is domi- 
nant only in the local surroundings of the MBH. Unlike stel- 
lar feedback, which operates globally, this new suppression 
mechanism does not require additional star formation and/or 
extensive mass expulsion out of the galactic potential. 

4.4. Regulated Black Hole Growth 

Heating by MBH feedback, which locally suppresses star 
formation, also makes the MBH to self-regulate its own 
growth. Figure Q3] shows the MBH accretion history versus 
time for Sim-SF and Sim-RMF. 

• The MBH in Sim-SF has grown exponentially to 3 x 
10 6 M Q in 350 Myr, already about >30 times more mas- 
sive than its initial mass of 10 5 M o . The MBH main- 
tained the accretion rate of 0.2 - 0.6 MEdd during this 
time period, corresponding to the unhindered growth of 
the MBH when there is no mechanism to self-regulate 
itself other than stellar feedback. 

• Over the course of the same period the MBH in Sim- 
RMF has grown by only ~ 70% to 1.7 x 1O 5 M . The 
heated and diffused ISM in the vicinity of the MBH 
considerably suppresses the Bondi-Hoyle accretion es- 
timate to as low as ~ 0.02 M^dd m this period. This 
indicates that the MBH feedback described in previous 



sections is a possible "radio-m ode" analogue where the 
accretion rate is ~ 0.05M E dd (ICroton et al.ll2006h . It 
also presents a potential route to the relatively low mass 
MBH at the center of the Milky Way (3 x 1O 6 M ). 20 

Therefore, the feedback from a MBH is confirmed as an 
effective mechanism for slowing down the accretion of gas 
onto itself. Without having to suddenly unbind all the sur- 
rounding gas, the MBH self-regulates its growth by heating 
up the neighborhood and keeping it hot for an extended pe- 
riod of time. This finding is consistent with the work bv De- 
buhr et al. d2011h . who claimed that the growth of a MBH 
could be "self-regula t ed", rather than "s upply-limited" (as in 
Sprin gel et afll2005bt iTeyssier et al.ll2O10h where quasar-like 
MBH feedback drive energetic large-scale outflows to unbind 
a significant amount of gas. 

5. DISCUSSION 

5.1. Summary and Conclusions 

A state-of-the-art numerical framework which fully incor- 
porates gas, stars, and a central massive black hole is devel- 
oped. Our simulation, for the first time, followed the compre- 

Authors again caution that the initial MBH mass (10 5 M ra ) and fina l 
masses in both Sim-SF and Sim-RMF lie below the Mag orrian etal] (1998) 
relationship, when 10% of the stellar mass is assumed to be in the bulge. 
The choice of a relatively small initial MBH mass may have caused a "radio- 
mode"-like MBH feedback, and a slower mass growth of the black hole. 
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FIG. 15. — Top: black hole mass accretion history. Note that the mass of 
the MBH of Sim-RMF has not doubled during this period, while the MBH of 
Sim-SF has grown exponentially. Bottom: mass accretion rate onto the MBH 
in the unit of Eddington rate. In each panel, the red solid line and the green 
dashed line represent Sim-SF and Sim-RMF, respectively. 

hensive evolution of a massive star-forming galaxy with self- 
consistently modeled stars and a MBH. Our novel framework 
renders a completely different, yet physically more accurate 
picture of how a galaxy and its embedded MBH evolve under 
each other's influence, providing a powerful means in under- 
standing the coevolution of galaxies and MBHs. Our main 
results and new advancements are as follows. 

1 . Molecular Cloud Formation and Feedback: We have 
included a new model of molecular cloud formation 
and stellar feedback in our code ( ^2.4l to l2~5l l. Unlike 
previous star formation recipes based on the Schmidt 
relation, a particle spawns when a gas cell of a typical 
molecular cloud size, 15.2 pc, actually becomes Jeans 
unstable. Then the molecular cloud particle gradually 
produces stellar mass while returning a large fraction 
of mass back to the gas with thermal feedback energy, 
modeling the observed slow star formation in molecu- 
lar clouds. Thermal stell ar fe edback is shown to self- 
regulates star formation ( 34.11 . 

2. Massive Black Hole Accretion and Feedback: We have 
successfully developed a self-consistent model of ac- 
cretion of gas onto a MB H and i ts ra diative and me- 
chanical feedback effects ( £|2.6l to |2~8T >. Gas accretion 
onto the MBH is estimated with the Bondi-Hoyle for- 
mula, but without any boost factor, as we begin to re- 
solve the Bondi radius. Monochromatic X-ray photons 
from the MBH are followed through three dimensional 
adaptive ray tracing, rendering the radiative feedback 



of a MBH; here, rays of photons ionize and heat the 
gas, and exert momentum onto the gas. Finally, the 
mechanical feedback of the MBH is represented by 
bipolar jets with velocities of ~ 10 3 kms _1 launched 
from the vicinity of the MBH accretion disk, well re- 
solved in our high-resolution AMR simulations. Our 
approach is significantly different from the previous 
recipes for MBH feedback in galactic scale simulations 
to date (e.g. ISprineel et al.ll2005bh ISiiacki et al] 120071: 
iBooth & Schavd 120091: iTevssier et al J 120101: Debuhr 
et al. 1201 ll) . vet more accurately presents the physics 
of MBHs when properly incorporated with a high dy- 
namic range. 

3. Locally Suppressed Star Formation: By investigating 
the coevolution of a 9.2 x 10 n MQ galactic halo and its 
10 5 M e embedded MBH at z ~ 3, we show that MBH 
feedback, especially its radiation, heats the surrounding 
ISM up to 10 6 K through photoionization and Compton 
heating and thus locally suppresses star formation in the 
inner core of a galaxy ( 34.21 i. The feedback also con- 
siderably changes the stellar distribution at the galactic 
center. This new channel of feedback from a slowly 
growing MBH is particularly interesting because it is 
only locally dominant, and does not require the heating 
of gas globally on the disk, or instigat e a m assive gas 
expulsion out of the galactic potential ( £|4.31 l. 

4. Self-regulated Black Hole Growth: MBH feedback is 
also demonstrated to be an effective mechanism for 
slowing down the accretion of gas onto the MBH it- 
self. Without necessarily unbinding all of its surround- 
ing gas, the MBH self-regulates its growth by keeping 
the surrounding ISM hot for an extended period of time 
( £|4.41 >. Therefore, our results possibly are consistent 
with a "radio-mode" analogue of MBH feedback. 

Our method limits the use of ad hoc formulation and instead 
more accurately models the physics of galaxy formation. As 
a result, four key components of galactic scale physics, (a) 
molecular cloud formation, (b) stellar feedback, (c) MBH ac- 
cretion, and (d) MBH feedback, work self-consistently in one 
comprehensive framework. As an example, the radiation and 
jets from the MBH heat up the surrounding gas and create 
hot regions, but the thermal couplings of the radiative and 
mechanical energy are all carried out by the shock-capturing 
radiation hydrodynamics AMR scheme itself, not by any pre- 
supposed thermal deposition model. In our framework, one 
should also be able to couple small-scale physics (such as 
molecular cloud formation and feedback) with large-scale 
physics (such as quasar-driven galactic outflows) without any 
sub-resolution model. These first results undoubtedly demon- 
strate that we can now develop an unabridged, self-consistent 
numerical framework for both galaxies and MBHs. 

5.2. Future Work 

While proven to be fruitful already in producing robust re- 
sults, our comprehensive galaxy formation framework is only 
the first step forward in the right direction. Imminent future 
projects and improvements are as follows. 

1 . Merging Galaxies: What is experimented in this work 
is a quiescent form of MBH feedback, possibly a 
"radio-mode" analogue. Meanwhile, MBH feedback is 
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known to make a dramatic difference in galaxy merg- 
ers, which is frequent in hierarchical structure forma- 
tion. The gas funneled into gala ctic centers trigger s 
"quasar-mode" MBH feedback dHopkins et all 12008). 
which subsequently reshapes the relation between black 
hole masses and bulge stellar distributions (Di Matteo 
et al. 12005b iJohansson et alj|2009l) . Detailed research 
on this more intense mode of MBH feedback will be 
the topic of a future paper. 

2. Parameter Studies: More comprehensive parameter 
studies should follow, especially in the parametrization 
of MBH feedback and the efficiency of stellar feedback. 
The results should be compared and calibrated with ob- 
servations such as bul ge to disk mass rati o and gas to 
stellar mass ratio (e.g. lGavazzi et al.ll2008h. or wi th an- 
alytical investigations (e.g. ISazonov et aU l2005h . In 
particular, the disk-bulge decomposition of simulated 
galaxies will be the subject of subsequent analysis of 
our simulations. 

3. Improving Mechanical Feedback: In the results pre- 
sented herein, mechanical feedback is energetically 
secondary to radiative feedback because the mass ac- 
creted onto the MBH is not large enough to repeatedly 
drive jets. These infrequent jets easily penetrate the 
ISM without necessarily creating sizable shocks or en- 
training a large amount of gas. However, a few mech- 
anisms will be considered in the future which could 
have enhanced the effectiveness of jets. Magnetic fields 
could aid the jets in efficiently depositing outflow mo- 
mentum onto the infalling gas, as was shown by the 
studies on t he evolution of jets in the presence o f mag- 
netic fields dDubois et al.ll2009l : IWang et alJ2010T) . Cos- 
mic rays accelerated by relativistic jets and shock fronts 
(Skillman et al. 2008) could boost the effectiveness of 
jets, too. 

4. Improving Radiative Feedback: For now, monochro- 
matic X-ray photons are utilized to carry the energy of 
MBH radiation ( fl2.7l : however, a better model will be 
needed to describe the polychromatic energy distribu- 
tion of MBH radiation. Ideally one wants to have a 
large number of spectral energy bins, each of which 
is separately followed through three-dimensional ray 
tracing. Given the computationally challenging nature 
of polychromatic radiative transfer, however, tabulated 
rates of photoionization and photoheating as functions 
of optical depth can be a good alternative. Moreover, to 
accurately quantify the radiative feedback on the gas in 
the vicinity of a MBH, the pressure force on dust grains 
needs to be computed. This could have increased the 
radiation pressure in the presented results, especially in 
th e central < kpc regi on. For this purpose, dust models 
in lRocha"e t al. (2008) will need to be considered. 

5. Adding Supplementary Feedback Channels: A MBH 
particle in our work represents not just the black hole 
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itself, but also includes accreting gas and stars deep 
within the galactic nucleus. Thus, there is a need for 
other feedbac k channels, such as stellar winds from a 
nuclear disk (ICiotti et al.l2009h . The nuclear disk winds 
can be implemented as thermal deposition of energy, 
working in conjunction with the aforementioned radia- 
tive and mechanical feedback. Stellar UV radiation 
from the nuclear disk can also be incorporated into the 
radiative feedback of the MBH. Including this supple- 
mentary feedback will reveal the multi-faceted nature 
of the coupling of MBH energy with its surroundings. 

6. Improving Accretion Estimate: The accretion estimate 
using the Bondi-Hoyle formula will need to be im- 
proved, especially when the gas disk around the MBH 
can be resolved down to the Bondi radius. Different 
estimates such as the ones considering gas angular mo- 
mentum (Hopkins & Ouataerj2010l : lLevine et alT2 010) 
are attractive candidates that should be explored. 

7. Nonthermal Pressure Sources: Nonthermal pressur e 
sources such as magnetic fields (Wa ng & Abelll2009h . 
stellar UV radiation, and cosmic rays are missing in this 
work, but should be included in future simulations. 

We also recognize that the results from our experiment can 
provide the community with better sub-resolution models for 
MBH physics. For example, the radial profile of heating 
rates by the MBH in our simulation can be tabulated; in a 
coarsely resolved particle-based simulation, one can deposit 
thermal energy according to this radial dependence into a vol- 
ume larger than a typical smoothing kernel. This can be a 
useful means for improving the particle-based simulations as 
well as for speeding up future large-scale AMR calculations, 
such as the formation of high-redshift quasars and the reion- 
ization of intergalactic helium. 
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APPENDIX 
A. ACCRETION RATE ESTIMATE 
The Bondi-Hoyle accretion rate estimate at the molecular cloud formation threshold n = n t hres = 125 cm~ 3 and c s = 10 kms -1 

is 

Mb — 0.004 (Mbh/1O 5 M0) 2 M©yr~', (Al) 

which is bigger than the Eddington rate, 

M Edd ~ 0.002 (M BH /10 5 M o ) M q yr" 1 . (A2) 

Therefore the density threshold for molecular cloud formation does not limit the accretion rate at any time. To put it in another 
way, because the Bondi-Hoyle rate can surpass the Eddington rate in dense clumps of n ~ «thres> the Eddington limit should play 
a crucial role in restricting the accretion. 

We note that, even without the aid of the boost factor (unlike in Eq.([TJ), the Bondi-Ho yle e stimate in our simulations can 
surpass the Eddington limit, and averages at 0.2 - 0.6 M^Ad m the reported simulation (See £|4.41 i. In other words, had we used 
the boost factor the Bondi-Hoyle estimate would have been almost always limited by the Eddington limit. It is partly because 
the gas density in our simulations reaches up to n = »thres m the finest cells. Th is justifies our choice of not employing the boost 
factor which has been common in the previous work (e.g. Spri ngel et al.ll2005bl and many others). However, the omission of the 
usual boost factor does not indicate that our calculation captures the turbulent accreting flow around the MBH accretion disk. No 
contemporary galactic scale simulation - including the reported simulation - has ever captured the turbulent interstellar medium 
which exists well below the typical resolution limit. Therefore, many other models for the MBH accretion estimate are equally 
applicable in galactic sc ale simulations, including the Bondi-Hoyle estimate with a boost factor parametrized by the gas density 
lBooth&Schavel(l2009l) . 

Readers should keep in mind that the spherical Bondi-Hoyle estimate or any of its variations provides only an approximated 
accretion rate in a low resolution simulation. In reality it is unlikely that the gas around the accreting black hole is spherically 
symmetric. Rather, the gas is expected to reside in a rotationally supported disk which is far from being resolved in a prese nt-day 
galactic scale simula tion. In this regard, attempts are being made to improve the Bond i-Hoyle estimate dLevine et al.ll2Q10f) or to 
develop alternatives dHopkins & Ouataerj|20ToLlPower et al.ll2010l:lDebuhr et aTll201 lh . 



